subroutine read_matrix(A, B)
  implicit none
  integer, parameter :: dbl = 8, fnum = 4
  integer :: dim, b_col
  real(dbl), allocatable, dimension(:,:), intent(inout) :: A
  real(dbl), allocatable, dimension(:,:), intent(inout) :: B
  integer :: k, status

  open(unit=fnum, file='matrix.txt', status='old', action='read', iostat=status)

  openif: if(status /= 0) then
     write(*,*) 'open file error'
     stop
  else
     read(fnum,*) dim, b_col
     allocate(A(dim,dim))
     allocate(B(dim,b_col))

     readA: do k = 1,dim
        read(fnum,*) A(k,:)
     end do readA

     readB: do k = 1, b_col
        read(fnum,*) B(:,k)
     end do readB

  end if openif

  close(fnum)
end subroutine read_matrix
